
import sys
import featureindex

def parseMatepairs(file1, file2):

    # record second-matepair locations
    matepairRanges = {}
    n = 0
    for feature in featureindex.getFeatures(file2):
        details = feature.getDetails()
        if not details.id.endswith("/2"):
            raise Exception("bad sequence id: ", id)
        id = details.id.replace("/2", "")
        spos = long(details.spos)
        epos = long(details.epos)
        chr = details.chr
        strand = details.strand
        seq = details.seq
        if matepairRanges.has_key(id):
            # this is the second half of a splice junction, so merge seq entries
            loc = matepairRanges.get(id)
            loc[2] = min(spos, loc[2])
            loc[3] = max(epos, loc[3])
            if strand=="+":
                loc[4] = seq + loc[4]
            else:
                loc[4] = loc[4] + seq
            matepairRanges[id] = loc
        else:
            matepairRanges[id] = [chr, strand, spos, epos, seq]
        n += 1
        if 0==(n%100000):
            print >> sys.stderr, "# Read ", n, " features from ", file2

    # walk through first-matepair seqs
    matepairOneRanges = {}
    n = 0
    for feature in featureindex.getFeatures(file1):
        details = feature.getDetails()
        if not details.id.endswith("/1"):
            raise Exception("bad sequence id: ", id)
        id = details.id.replace("/1", "")
        spos = long(details.spos)
        epos = long(details.epos)
        chr = details.chr
        strand = details.strand
        seq = details.seq
        if details.startDrop != "0" or details.endDrop != "0":
            # first-matepair seq mapped to splice junction
            if matepairOneRanges.has_key(id):
                # this is the second half, so handle it now
                prevLoc = matepairOneRanges.get(id)
                spos = min(spos, prevLoc[2])
                epos = max(epos, prevLoc[3])
                # combine sequence fragments from splice junction:
                if strand=="+":
                    seq = prevLoc[4] + seq
                else:
                    seq = seq + prevLoc[4]
            else:
                # store this half and wait for second half to be seen
                matepairOneRanges[id] = [chr, strand, spos, epos, seq]
                continue
        oldLoc = matepairRanges.get(id, None)
        if oldLoc == None:
            print "%s\tFIRST\t%s\t%s\t%ld\t%ld\t%s\t%s\t%s\t%s\t%s\t%s" % (id, chr, strand, spos, epos, seq, "NA", "NA", "NA", "NA", "NA")
        else:
            if oldLoc[0]==chr and oldLoc[1]==strand:
                t = "PAIR"
            elif oldLoc[0]==chr and oldLoc[1]!=strand:
                t = "PAIR-STRAND"
            else:
                t = "PAIR-CHR"
            print "%s\t%s\t%s\t%s\t%ld\t%ld\t%s\t%s\t%s\t%ld\t%ld\t%s" % (id, t, chr, strand, spos, epos, seq, oldLoc[0], oldLoc[1], oldLoc[2], oldLoc[3], oldLoc[4])
            del matepairRanges[id]
        n += 1
        if 0==(n%100000):
            print >> sys.stderr, "# Read ", n, " features from ", file1

    # handle seqs that mapped to second matepair only
    for id in matepairRanges.keys():
        oldLoc = matepairRanges.get(id)
        print "%s\tSECOND\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%ld\t%ld\t%s" % (id, "NA", "NA", "NA", "NA", "NA", oldLoc[0], oldLoc[1], oldLoc[2], oldLoc[3], oldLoc[4])

if __name__=="__main__":
    parseMatepairs(sys.argv[1], sys.argv[2])
